#### 7.2: Consistency with Existing Findings on US Presidential Hawkishness ####

source("./code/loadPackages.R")

#### Replicating Horowitz and Stam (2014) ####

## Horowitz and Stam data
# MIDs data
hs = haven::read_dta("./otherdata/HorowitzStamLeadersIOMIDReplication.dta")
hsu = hs[which(hs$ccode==2),]

# War data
hsw = haven::read_dta("./otherdata/HorowitzStamLeadersIOWarReplication.dta")
hswu = hsw[which(hsw$ccode==2),]

# Simple model of MIDs initiation
hsFit0 = glm(cwinit ~ milnoncombat + combat, data=hsu, family="binomial")
hsFit0r = coeftest(hsFit0, vcov = vcovCL, cluster = ~ leaderid)
hsFit0r

# Full model of MIDs initiation
hsFit1 = glm(cwinit ~ milnoncombat + combat + rebel + warwin + warloss + rebelwin + rebelloss + 
               age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag + 
               cwpceyrs1 + cwpceyrs2 + cwpceyrs3, data=hsu, family="binomial")
hsFit1r = coeftest(hsFit1, vcov = vcovCL, cluster = ~ leaderid)
hsFit1r

# Simple model of war initiation
hsFit2 = glm(initiation ~ milnoncombat + combat, data=hswu, family="binomial")
hsFit2r = coeftest(hsFit2, vcov = vcovCL, cluster = ~ leaderid)
hsFit2r

# Full model of war initiation
hsFit3 = glm(initiation ~ milnoncombat + combat + rebel + warwin + warloss + rebelwin + rebelloss + 
               age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag + 
               cwpceyrs1 + cwpceyrs2 + cwpceyrs3, data=hswu, family="binomial")
hsFit3r = coeftest(hsFit3, vcov = vcovCL, cluster = ~ leaderid)
hsFit3r

#### Table A34: Replication of Models in Table 1 of Horowitz and Stam (2014), Only US Presidents ####
# (Note: Code immediately below produces underlying table; code under that produces results with clustered standard errors)
stargazer(hsFit0, hsFit1, hsFit2, hsFit3, column.sep.width = "-10pt", no.space=T, align=T, 
          omit.stat = c("ll", "rsq", "adj.rsq", "aic", "f", "ser"),
          omit = c("cwpceyrs"),
          covariate.labels = c("Military service, no combat",
                               "Military service, combat",
                               "Rebel service",
                               "Prior war win",
                               "Prior war loss", 
                               "Prior rebel win",
                               "Prior rebel loss",
                               "Leader age",
                               "Autocracy",
                               "Material capabilities",
                               "Tau B with system leader",
                               "Time in office",
                               "Five-year MID challenge lag"))
stargazer(hsFit0r, hsFit1r, hsFit2r, hsFit3r, column.sep.width = "-10pt", no.space=T, align=T,
          omit = c("cwpceyrs"),
          covariate.labels = c("Military service, no combat",
                               "Military service, combat",
                               "Leader age",
                               "Material capabilities",
                               "Tau B with system leader",
                               "Time in office",
                               "Five-year MID challenge lag"))


#### Replicating Carter and Smith (2020) ####

# Load and install some more packages
pacman::p_load(tidyverse,haven,rstan,shinystan,rstanarm,loo,
               janitor,sjlabelled,tidylog,tidybayes,apsrtable,
               pscl,haven,broom)

## Loading Data Sets ##
LEAD_Dataset_Article <- read_dta("./otherdata/CarterSmith/LEAD_Dataset_Article.dta")
LEADYear             <- read_dta("./otherdata/CarterSmith/WhyLeadersFightMonadicReplication_updated.dta")
DebsGoemans          <- read_dta("./otherdata/CarterSmith/DebsGoemans_leaders.2.9.dta")
AllThetas            <- read.csv("./otherdata/CarterSmith/AllThetas.csv")

## Adapting the main LEAD data ##
LEAD_Dataset_Article  %>%
  as_tibble()         %>% 
  remove_all_labels() %>%
  clean_names()       %>% 
  rename(male=gender) %>% 
  mutate(entry = na_if(entry, -666),
         irregular=case_when(entry == 1 ~ 1, 
                             entry == 0 ~ 0, 
                             TRUE ~ NA_real_),
         scienceengineer=case_when(science == 1 | engineering == 1 ~ 1, 
                                   science == 0 & engineering == 0 ~ 0, 
                                   TRUE ~ NA_real_),
         creative=case_when(writer == 1 | filmmusic == 1 ~ 1, 
                            writer == 0 & filmmusic == 0 ~ 0, 
                            TRUE ~ NA_real_),
         deathyear = na_if(deathyear, -777),
         xrreg     = na_if(xrreg,     -77 | -88),
         xropen    = na_if(xropen,    -77 | -88),
         xconst    = na_if(xconst,    -77 | -88),
         parreg    = na_if(parreg,    -77 | -88),
         parcomp   = na_if(parcomp,   -77 | -88),
         democracy = na_if(democracy, -77 | -88),
         autocracy = na_if(autocracy, -77 | -88),
         polity    = na_if(polity,    -77 | -88),
         exit      = na_if(exit,      -888)) -> LeadDataTemp

DebsGoemans %>% rename(leadid29=leadid) %>% 
  dplyr::select(leadid29,initiator) -> DG

LeadDataTemp %>% dplyr::select(leaderid,leadid29,milservice) -> thetasmerge
thetas <- as_tibble(cbind(AllThetas,thetasmerge))

LEADYearAllThetas <- inner_join(as_tibble(LEADYear), thetas)

LEADYearAllThetas  %>% filter(ccode==2) %>% 
  dplyr::select(leadername, leadname=leaderarchigos, year, cwinit,force2dv,milservice,theta1,theta2,theta3,theta4) %>% 
  drop_na() -> middatatests

DebsGoemansAllThetas <- inner_join(DG,thetas)

DebsGoemansAllThetas %>% 
  dplyr::select(leadid29, initiator,milservice,theta1,theta2,theta3,theta4) %>% 
  drop_na() -> icbdatatests

## Main US presidents we care about ##
uskey = unique(DebsGoemans$leadid[which(DebsGoemans$ccode==2)])
icbdatatests = icbdatatests %>% filter(leadid29 %in% uskey)


## Redo Carter and Smith's MID analysis with Horowitz and Stam's controls added
csu = middatatests
csu = csu %>% arrange(year)
csu$leaderyear = paste(csu$leadname, csu$year, sep=" ")

hsus = hsu %>% arrange(year, leaderid)
hsus$leaderyear = paste(hsus$leaderarchigos, csu$year, sep=" ")
hswu = hswu %>% arrange(year, leaderid)
hswu$leaderyear = paste(hswu$leaderarchigos, hswu$year, sep=" ")

hsucontrol = hsus %>% dplyr::select(rebel, warwin, warloss, rebelwin, rebelloss, age, 
                                    aut, cinc, tau_lead, officetenure1000, fiveyearchallengelag, 
                                    cwpceyrs1, cwpceyrs2, cwpceyrs3)

# Final data for Carter and Smith MIDs
csud = cbind(csu, hsucontrol)

# Final data for Horowitz and Stam wars including Carter and Smith hawk measures
csuformerge = csu %>% dplyr::select(leaderyear, theta1:theta4) %>% distinct()
hsuwc = merge(hswu, csuformerge, by="leaderyear", all.x=T)


# Redo MID initiation analysis
cwinitbase1c <- glm(cwinit ~ milservice + rebel + warwin + warloss + rebelwin + rebelloss + 
                      age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)
cwinittheta1c <- glm(cwinit ~ theta1 + rebel + warwin + warloss + rebelwin + rebelloss + 
                       age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)
cwinittheta2c <- glm(cwinit ~ theta2 + rebel + warwin + warloss + rebelwin + rebelloss + 
                       age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)
cwinittheta3c <- glm(cwinit ~ theta3 + rebel + warwin + warloss + rebelwin + rebelloss + 
                       age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)
cwinittheta4c <- glm(cwinit ~ theta4 + rebel + warwin + warloss + rebelwin + rebelloss + 
                       age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)

#### Table A35: Extension of Carter and Smith (2020): MID Initiation, Only US Presidents ####
stargazer(cwinitbase1c, cwinittheta1c, cwinittheta2c, cwinittheta3c, cwinittheta4c, 
          column.sep.width = "2pt", no.space=T, align=T, dep.var.labels = "MID initiation",
          omit = c("cwpceyrs"),
          covariate.labels = c("Military", "M1", "M2", "M3", "M4",
                               "Rebel service",
                               "Prior war win",
                               "Prior war loss", 
                               "Prior rebel win",
                               "Prior rebel loss",
                               "Leader age",
                               "Autocracy",
                               "Material capabilities",
                               "Tau B with system leader",
                               "Time in office",
                               "Five-year MID challenge lag"))

# Redo severe MID initiation
force2dvbase1c <- glm(force2dv ~ milservice + rebel + warwin + warloss + rebelwin + rebelloss + 
                        age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)
force2dvtheta1c <- glm(force2dv ~ theta1 + rebel + warwin + warloss + rebelwin + rebelloss + 
                         age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)
force2dvtheta2c <- glm(force2dv ~ theta2 + rebel + warwin + warloss + rebelwin + rebelloss + 
                         age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)
force2dvtheta3c <- glm(force2dv ~ theta3 + rebel + warwin + warloss + rebelwin + rebelloss + 
                         age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)
force2dvtheta4c <- glm(force2dv ~ theta4 + rebel + warwin + warloss + rebelwin + rebelloss + 
                         age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, family=binomial(link="logit"), data=csud)

#### Table A36: Extension of Carter and Smith (2020): Severe MID Initiation, Only US Presidents ####
stargazer(force2dvbase1c, force2dvtheta1c, force2dvtheta2c, force2dvtheta3c, force2dvtheta4c, 
          column.sep.width = "2pt", no.space=T, align=T, dep.var.labels = "Severe MID initiation",
          covariate.labels = c("Military", "M1", "M2", "M3", "M4",
                               "Rebel service",
                               "Prior war win",
                               "Prior war loss", 
                               "Prior rebel win",
                               "Prior rebel loss",
                               "Leader age",
                               "Autocracy",
                               "Material capabilities",
                               "Tau B with system leader",
                               "Time in office",
                               "Five-year MID challenge lag"))

## Redo Horowitz and Stam using Carter and Smith measures
hsFit3z = glm(initiation ~ I(milnoncombat|combat) + rebel + warwin + warloss + rebelwin + rebelloss + 
                age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, data=hsuwc, family="binomial")
summary(hsFit3z)
hsFit3zr = coeftest(hsFit3z, vcov = vcovCL, cluster = ~ leaderid)
hsFit3zr

hsFit3a = glm(initiation ~ theta1 + rebel + warwin + warloss + rebelwin + rebelloss + 
                age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, data=hsuwc, family="binomial")
summary(hsFit3a)
hsFit3ar = coeftest(hsFit3a, vcov = vcovCL, cluster = ~ leaderid)
hsFit3ar

hsFit3b = glm(initiation ~ theta2 + rebel + warwin + warloss + rebelwin + rebelloss + 
                age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, data=hsuwc, family="binomial")
summary(hsFit3b)
hsFit3br = coeftest(hsFit3b, vcov = vcovCL, cluster = ~ leaderid)
hsFit3br

hsFit3c = glm(initiation ~ theta3 + rebel + warwin + warloss + rebelwin + rebelloss + 
                age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, data=hsuwc, family="binomial")
summary(hsFit3c)
hsFit3cr = coeftest(hsFit3c, vcov = vcovCL, cluster = ~ leaderid)
hsFit3cr

hsFit3d = glm(initiation ~ theta4 + rebel + warwin + warloss + rebelwin + rebelloss + 
                age + aut + cinc + tau_lead + officetenure1000 + fiveyearchallengelag, data=hsuwc, family="binomial")
summary(hsFit3d)
hsFit3dr = coeftest(hsFit3d, vcov = vcovCL, cluster = ~ leaderid)
hsFit3dr

#### Table A37: Extension of Carter and Smith (2020): Interstate War, Only US Presidents ####
# (Note: Code immediately below produces underlying table; code under that produces results with clustered standard errors)
stargazer(hsFit3z, hsFit3a, hsFit3b, hsFit3c, hsFit3d, column.sep.width = "-10pt", no.space=T, align=T, 
          covariate.labels = c("Military", 
                               "M1", "M2", "M3", "M4",
                               "Rebel service",
                               "Prior war win",
                               "Prior war loss", 
                               "Prior rebel win",
                               "Prior rebel loss",
                               "Leader age",
                               "Autocracy",
                               "Material capabilities",
                               "Tau B with system leader",
                               "Time in office",
                               "Five-year MID challenge lag"))
stargazer(hsFit3zr, hsFit3ar, hsFit3br, hsFit3cr, hsFit3dr, column.sep.width = "-10pt", no.space=T, align=T,
          covariate.labels = c("Military",
                               "M1", "M2", "M3", "M4",
                               "Prior war win",
                               "Leader age",
                               "Material capabilities",
                               "Tau B with system leader",
                               "Time in office",
                               "Five-year MID challenge lag"))

#### Table A38: Replication of Table 2 in Carter and Smith (2020) Using Only Observations from the United States ####

## Panel A: ICB initiation
initiatorbase1 <- glm(initiator ~ milservice,family=binomial(link="logit"), data=icbdatatests)
initiatortheta1 <- glm(initiator ~ theta1,family=binomial(link="logit"), data=icbdatatests)
initiatortheta2 <- glm(initiator ~ theta2,family=binomial(link="logit"), data=icbdatatests)
initiatortheta3 <- glm(initiator ~ theta3,family=binomial(link="logit"), data=icbdatatests)
initiatortheta4 <- glm(initiator ~ theta4,family=binomial(link="logit"), data=icbdatatests)

# Column 1 (Military vs. other measures)
vuong(initiatorbase1,initiatortheta1)
vuong(initiatorbase1,initiatortheta2)
vuong(initiatorbase1,initiatortheta3)
vuong(initiatorbase1,initiatortheta4)

# Column 2 (M1 vs. other measures)
vuong(initiatortheta1,initiatortheta2)
vuong(initiatortheta1,initiatortheta3)
vuong(initiatortheta1,initiatortheta4)

# Column 3 (M2 vs. other measures)
vuong(initiatortheta2,initiatortheta3)
vuong(initiatortheta2,initiatortheta4)

# Column 4 (M3 vs. other measures)
vuong(initiatortheta3,initiatortheta4)


## Panel B: MID initiation
cwinitbase1 <- glm(cwinit ~ milservice,family=binomial(link="logit"), data=middatatests)
cwinittheta1 <- glm(cwinit ~ theta1,family=binomial(link="logit"), data=middatatests)
cwinittheta2 <- glm(cwinit ~ theta2,family=binomial(link="logit"), data=middatatests)
cwinittheta3 <- glm(cwinit ~ theta3,family=binomial(link="logit"), data=middatatests)
cwinittheta4 <- glm(cwinit ~ theta4,family=binomial(link="logit"), data=middatatests)

# Column 1 (Military vs. other measures)
vuong(cwinitbase1,cwinittheta1)
vuong(cwinitbase1,cwinittheta2)
vuong(cwinitbase1,cwinittheta3)
vuong(cwinitbase1,cwinittheta4)

# Column 2 (M1 vs. other measures)
vuong(cwinittheta1,cwinittheta2)
vuong(cwinittheta1,cwinittheta3)
vuong(cwinittheta1,cwinittheta4)

# Column 3 (M2 vs. other measures)
vuong(cwinittheta2,cwinittheta3)
vuong(cwinittheta2,cwinittheta4)

# Column 4 (M3 vs. other measures)
vuong(cwinittheta3,cwinittheta4)


## Panel C: Severe MID initiation
force2dvbase1 <- glm(force2dv ~ milservice,family=binomial(link="logit"), data=middatatests)
force2dvtheta1 <- glm(force2dv ~ theta1,family=binomial(link="logit"), data=middatatests)
force2dvtheta2 <- glm(force2dv ~ theta2,family=binomial(link="logit"), data=middatatests)
force2dvtheta3 <- glm(force2dv ~ theta3,family=binomial(link="logit"), data=middatatests)
force2dvtheta4 <- glm(force2dv ~ theta4,family=binomial(link="logit"), data=middatatests)

stargazer(force2dvbase1, force2dvtheta1, force2dvtheta2, force2dvtheta3, force2dvtheta4, 
          column.sep.width = "2pt", no.space=T, align=T, dep.var.labels = "Severe MID initiation",
          covariate.labels = c("Military", "M1", "M2", "M3", "M4"))

# Column 1 (Military vs. other measures)
vuong(force2dvbase1,force2dvtheta1)
vuong(force2dvbase1,force2dvtheta2)
vuong(force2dvbase1,force2dvtheta3)
vuong(force2dvbase1,force2dvtheta4)

# Column 2 (M1 vs. other measures)
vuong(force2dvtheta1,force2dvtheta2)
vuong(force2dvtheta1,force2dvtheta3)
vuong(force2dvtheta1,force2dvtheta4)

# Column 3 (M2 vs. other measures)
vuong(force2dvtheta2,force2dvtheta3)
vuong(force2dvtheta2,force2dvtheta4)

# Column 4 (M3 vs. other measures)
vuong(force2dvtheta3,force2dvtheta4)